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Flow in a pressurized, vapor condensation (PVC) boron nitride nanotube (BNNT) 
production rig is modeled. A laser provides a thermal energy source to the tip of a 
boron fiber bundle in a high pressure nitrogen chamber causing a plume of boron- 
rich gas to rise. The buoyancy driven flow is modeled as a mixture of thermally 
perfect gases (B, B 2 , N, N 2 , BN) in either thermochemical equilibrium or chemical 
nonequilibrium assuming steady-state melt and vaporization from a 1 mm radius 
spot at the axis of an axisymmetric chamber. The simulation is intended to de- 
fine the macroscopic thermochemical environment from which boron-rich species, 
including nanotubes, condense out of the plume. Simulations indicate a high tem- 
perature environment {T > 4400K) for elevated pressures within 1 mm of the surface 
sufficient to dissociate molecular nitrogen and form BN at the base of the plume. 
Modifications to Program LAURA, a finite- volume based solver for hypersonic flows 
including coupled radiation and ablation, are described to enable this simulation. 
Simulations indicate that high pressure synthesis conditions enable formation of 
BN vapor in the plume that may serve to enhance formation of exceptionally long 
nanotubes in the PVC process. 


I. Introduction 

The future of manned and unmanned space flight and exploration depends on economical access 
to space. A key enabler of economic access to space is the use of multifunctional lightweight 
materials. Boron nitride nanotube (BNNT) composites offer distinct advantages for enhanced 
survivability during long-term exploration flights. BNNTs are structurally a close analog of carbon 
nanotubes (CNTs), with the carbon atoms alternately substituted with boron and nitrogen atoms. 
BNNTs are durable at low temperature in space (-157 C) and exhibit good thermal stability up 
to 800 C with a low coefficient of thermal expansion. Their strength-to-weight ratios and stiffness 
are greater than 95% of carbon nanotube values with an elastic modulus approaching 1 TPa. (The 
fibers are calculated to be the world’s strongest structural fiber above 400 C.) Boron has one 
of the largest neutron capture cross sections of all elements in the periodic table. Nitrogen has 
a larger neutron capture cross section than carbon. Consequently, BNNT materials have great 
potential for radiation shielding applications. Their high strength also makes BNNT fabrics well 
suited for deployable aerobrake surfaces, offering dual use opportunities for radiation shielding in 
transit and deployable aerobrake on entry into a planetary atmosphere. These materials can enable 
development of aerospace structures with both impact resistance and radiation shielding capability, 
thereby improving reliability and durability for long-term missions. 
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Arenal et al. 1 have recently reviewed BNNT production methods and associated nanotube 
properties. Not included in this review is a description of a recently developed and unique high- 
temperature, pressurized vapor condensation (PVC) method. 2 The NASA Langley Research Center 
(LaRC), National Institute of Aerospace (NIA) and Jefferson Sciences Association/DoE team dis- 
covered that this novel synthesis method produces thin, long, high aspect ratio, highly crystalline 
BNNTs using a conventional laser at high pressure. 2 The PVC method appears to minimize tube 
defects and create exceptionally long (order mm) tubes. The synthesis of longer tubes is advan- 
tageous for structural reinforcement applications, where it is more likely that the desired strength 
properties will extend from the nanoscale to the macroscale. 

Macroscopic details of the PVC production process (plume formation, species densities, tem- 
peratures, flow rates) ensuing from laser energy deposition on a boron sample in a pressurized 
test rig are modeled herein. The simulations here initially assume axisymmetric, single phase flow. 
Neither the possible presence of boron droplets sputtered from the boron source nor the nucleation 
of solids from the plume are included in these initial simulations. The intent here is to guide a 
subsequent experimental program to measure the thermochemical environment in the plume. The 
more critical components of nanotube formation on the molecular dynamics level (homogeneous 
nucleation, heterogeneous nucleation, surface catalysis) will require input from the macroscopic 
environment to be understood. Both the macroscopic environment and the molecular dynamics 
must be accurately modeled in order to optimize BNNT production rates and ultimately create 
an industrial scale processing capability. Our objective is a coordinated modeling effort to address 
these challenges to enable an industrial level BNNT production. To this end, we focus on the flow 
rates of BN vapor as a function of pressure to possibly explain improved qualities of nanotubes 
produced with the PVC method. 


II. Nomenclature 



Roman symbols 


Greek symbols 

A r , B r 

Arrhenius rate parameters, Eq. 1 

a 

absorbed radiation fraction 

c 

mass fraction 

e 

emissivity 

9 

gravitational acceleration 

P 

density 

h 

enthalpy per unit mass 

a 

Stefan-Boltzmann constant 

A H f 

heat of formation 

X 

mole fraction 


heat of vaporization for boron 



k 

reaction rate coefficient 


Subscripts 

K c 

equilibrium constant 

0 

conditions in chamber 

m 

mass loss rate, mg/mm 2 -hr 

act 

activation 

fh 

integrated tin, mg/hr 

b 

backward rate 

n 

number of moles 

cond 

conduction 

P 

pressure 

conv 

convection 

q 

heating rate 

f 

forward rate 

r 

radius 

melt 

property of droplet 

R 

universal gas constant 

r 

reaction number 

t 

time 

s 

species 

T 

temperature 

sat 

saturation 

v n 

velocity normal to wall 

vap 

vapor 

z 

height above droplet 

w 

wall 
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III. Pressurized Vapor Condensation (PVC) Method - Process Description 

The PVC method has been described by Smith et al. 2,3 The process is summarized here so 
that the models and simulations may be better understood in relation to actual conditions. 

A. System Components 

A schematic of the PVC system for producing Boron Nitride NanoTubes (BNNT) is presented in 
Figure 1. A laser beam (red line) enters a pressurized chamber (component H) through a window 
(component G) at the top of the chamber. A feedstock of bundled boron fibers (component K) 
is slowly guided into the chamber using the target feed drive (component J) and target adjusters 
(component L). A nitrogen supply line (component I), flow controls (component O), and exhaust 
flow meter (component D) are all used to maintain a constant pressure environment in the chamber. 
The laser is focused toward the tip of the boron fiber bundle in the chamber. 



C. 


U U 



Figure 1. (Provided by Michael Smith, BNNT LLC) Schematic of the BNNT production rig using 
PVC method: A. Exhaust to building exterior; B. HEPA filter (0.3 micron); C. Sliding, bullet-proof, 
protective shell; D. Chamber exhaust flowmeter; E. Beam-shaping optics; F. Turning mirrors; G. ZnSe 
Input window; H. Pressure chamber; I. Nitrogen supply line; J. Target feed drive; K. Target (B or 
BN); L. Vertical/Horizontal target adjusters; M. Pico- motor-controlled turning mirror; N. C02 laser 
beam; O. Pressure flow controls; P. Electronic controls (target and beam manipulation); Q. Optics 
table; R. Condenser Loop; S. View port (video cameras). 


B. Formation of Molten Boron Droplet 

An illustration of the interaction of the laser and the boron fiber bundle (not to scale) is presented 
in Fig. 2. The laser power and focus is adjusted to strike the boron fiber bundle in a way that 
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enables a boron droplet to suspend by surface tension from the end of the fiber bundle. The laser 
then directly strikes only the droplet. Energy is conducted and re-radiated from the droplet into 
the fiber bundle. In a perfectly tuned system all mass loss is through vaporization, the melt rate to 
replenish the droplet exactly matches the mass loss rate through vaporization, and the fiber bundle 
is continuously fed into the chamber to keep the interaction location fixed within the chamber. 
If delivered power is too great the droplet boils and drops off the end of the fiber bundle. Any 
such perfectly tuned combination of delivered laser power and constant feed rate of the boron fiber 
bundle into the chamber without any droplets falling off is considered a steady state vaporization 
condition. It requires a balance of mass conversion from solid phase to liquid phase to gas phase. 
It requires a corresponding balance of energy to drive these phase changes from the laser to the 
droplet to the fiber bundle along with convective and radiative losses to the surrounding system. 
The intent of the simulations to be described in the next section is to model the steady state 
vaporization. 



Figure 2. Laser aimed at the tip of a boron fiber bundle causes molten droplet to suspend from the 
tip, held there by surface tension, and a plume of boron-rich gas rises from the hot spot driven by 
buoyancy. Boron nitride nanotubes form on a wire condenser loop (not shown) when the rising plume 
contacts it. 


C. Formation of Boron Nitride Nanotubes on Condenser Wire 

Also shown in Fig. 2 is a plume of hot gas fed by vaporization off the droplet. The plume rises by 
buoyancy in which the low density, hot gas in the plume base is displaced vertically by the higher 
density, cooler nitrogen gas surrounding it. As the nitrogen is drawn in to the base of the plume 
it is heated and reacts with the boron vapor. Not shown in Fig. 2 is a condenser wire loop which 
serves as a nucleation site for BNNTs and other boron species when the plume flows over it. The 
condenser wire loop (component R) in Fig. 1 can be pulled like a clothesline on pulleys within the 
chamber. This action draws condensate out of the way of the plume and exposes fresh nucleation 
sites. In the configuration shown in Fig. 1 the condensate is harvested after the system is powered 
off and the chamber can be de-pressurized and opened. Simulations of the gas phase in the plume 
are intended to define the constituents and thermodynamic environment available for nucleation. 
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IV. Modeling 


A. Computational Domain, Grid, and Boundary Conditions 

The discretization of an idealized chamber is presented in Fig. 3. The model assumes that the actual 
plume is axi-symmetric and that its evolution within the actual chamber up to the condenser wire 
is not strongly effected by the orientation of chamber walls far removed (50 to 100 plume diameters) 
from the plume axis. This assumption enables creation of an idealized, axi-symmetric chamber to 
simplify the simulation and enable parametric investigations in reasonable time. The axis of the 
chamber corresponds to the axis of the plume as seen in the center of the figure. The idealized 
chamber radius is 5 cm. The idealized chamber height is 10 cm. The radius of the laser spot in 
the droplet is set to 1 mm for all simulations. At time t = 0 the chamber is filled with molecular 
nitrogen at temperature To, pressure po, and density po with zero velocity throughout. 

The top and bottom boundaries, except for the 1 mm radius spot heated by the laser, is 
treated as a cold ( T w = 300. K), no-slip wall. Species mass fractions at the walls are computed 
assuming chemical equilibrium at the wall temperature and the locally computed elemental mass 
fractions from the elemental continuity equations. The left boundary is the axis of symmetry. The 
right boundary is treated as a porous wall to maintain constant pressure in the chamber in an 
environment in which mass and energy are added continuously. In the actual chamber nitrogen 
feed and purge lines (grey lines in Fig. 1) purge gas from the chamber and replace it with molecular 
nitrogen to maintain constant pressure. The porous wall boundary condition replaces this function 
while maintaining an axisymmetric flow environment. It computes the difference between the 
local value of pressure at the wall and the target pressure, A p w = p w — pt ar ■ It is derived from 
Bernoulli’s Equation and assumes an idealized porosity with zero losses. If A p. w > 0 then flow at 
local conditions is directed out of the chamber with velocity V n = — i/2A p w / p w . If A p w < 0 then 
pure nitrogen flow at To = 300 K is directed into the chamber with velocity V n = \J — 2Ap w /po. 

In the lower right inset of Fig. 3 the domain above the 1 mm radius irradiated spot is featured. In 
this simulation the droplet surface has no curvature. In fact, the largest and potentially significant 
difference between this idealized domain and the actual domain is that the idealized domain treats 
the irradiated surface of the droplet as a spot on the floor of the chamber. In the actual chamber, 
gas can rise from beneath the droplet suspended from the end of the fiber bundle. In the present 
simulations, gas flows along a cold wall before feeding the base of the plume. Possible consequences 
of this simplification are discussed in a later section. The cell size normal to the wall in the nominal 
grid is 10~ 7 nr. The tight spacing here helps maintain stability in the very early evolution of the 
plume when the spot is first irradiated. The cell size tangent to the wall over the 1 mm radius is 
a constant 2 1CP 5 nr. Mesh size grows in the radial direction by a factor of 1.1 beyond the 1 mm 
edge of the droplet. Mesh size grows by the same 1.1 factor in the vertical direction. Mesh growth 
is terminated when the cell dimension achieves 2 mm. This grid generation algorithm maintains a 
fine spacing in the radial direction in the core of the chamber to capture diffusion across the edge 
of the plume. Additional discussion of grid resolution will be presented in a later section on Grid 
Convergence. The boundary condition across the 1 mm radius at the bottom of the domain defines 
the surface response model. The formulation of this boundary condition is discussed subsequently 
in the section titled Material Response Model. 

B. Gas Properties 

The nominal atmosphere inside the pressurized chamber is 100% N 2 . Parameter studies also con- 
sider the effects of noble gas (He or Ar) in the chamber. The B-N system is assumed to include 
(B, BN, N, N 2 , and B 2 ) - all in the gas phase. An option for multiple phases in the plume is not 
yet considered and so the simulations are intended to describe conditions up to the first nucleation 
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Steady State 
Ablation 


Figure 3. Boundary conditions and grid (113 x 144 cells) used to discretize idealized, axisymmetric 
pressure chamber (5 cm radius x 10 cm height). Upper right inset shows perspective of 5 degree 
wedge shaped domain. Lower right inset shows detail above 1 mm radius spot where the laser heats 
the droplet at the tip of the boron fiber bundle. 


event. Curve fits for thermodynamic properties of these species are obtained from the CEA code. 4 
Transport properties are calculated as functions of collision cross-sections and temperature.' 5 Col- 
lision cross sections for the B-N system were approximated by using the corresponding collision 
cross sections in the C-N system 6-8 along with the appropriate molecular weights of the colliding 
pairs in the B-N system. A spot check with an evaluation of B-B cross sections 9 at 2000 K indicate 
they are approximately 20% higher than the corresponding C-C cross sections. 10 

In the case of flow in chemical equilibrium the elemental continuity equations are solved and a 
free energy minimization is implemented in every control volume to solve for species mass fractions 
as a function of elemental mass fraction, temperature, and pressure. 11 

Early equilibrium simulations had been run with only a four species model that neglected B 2 . 
These simulations indicated that the presence of B 2 significantly diminished the amount of BN in 
the upper part of the plume. A quick check on the possible influence of reaction path to equilibrium 
was then initiated. The intent here was to see how BN, created in the hot base of the plume, was 
consumed as the gas rises and cools. Is it credible that BN levels could freeze (remain relatively 
constant) as the plume rises and cools, even at elevated pressures? This question regarding the 
qualitative effects of kinetics in a rising plume is addressed by approximating a kinetic model for the 
B-N system with a C-N system template as shown in Table 1. (An experimentally verified kinetic 
model for the B-N system has not yet been identified.) Column 1 is the reaction number, Column 
2 shows a proposed reaction in the B-N system and column 3 shows the corresponding template 
reaction from the C-N system. Column 4 presents the pre-exponential factor of the Arhenius 
rate formulation as defined for the template reaction, 12,13 where the Arrhenius rate coefficients for 
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forward and backward reactions are given by 


kf, r = A r T Br exp(-T r , act /T) (1) 

kb,r = kf tr /K Ctr (T) (2) 

The activation temperature T r ^ ac t for reaction r is approximated using the heats of formation of 
reactants and products compiled in Table 2. 

( products reactants \ 

Y n s , r AH f , s - Y n a , r AH f A/R (3) 

Column 5 of Table 1 presents the proposed activation temperature for the B-N system reaction as 
calculated from Eq. 3. Column 6 presents the activation temperature used in the template C-N 
system and column 7 presents the value of T* as a reference condition for the credibility of the 
approximation in Eq. 3. T* is within 3% of the experimentally determined activation temperature 
for reactions 1-3. The approximation is not so good with reactions 4 and 5 where differences are 47% 
and -25%, respectively. The dissociation of N 2 by collision with any generic partner is included in 
the model. It does not appear in Table 1 because its Arrhenius rate coefficients are readily available 
and do not require a template. 


Table 1. Chemical Kinetic Models for B-N and C-N Systems 


r 

B-N system 
proposed reactions 

C-N system 
template reactions 

pre-expon. 

A r 

proposed 

T r ,act, K 

template 
T r ,act, K 

template 
T* , K 

1 

B 2 + M <=> 2B + M 

C 2 + M <=> 2C + M 

3.7 10 14 

34914 

69900 

71637 

2 

BN + M <=> B + N + M 

CN + M <=> C + N + M 

2.5 10 14 

66834 

87740 

90717 

3 

N 2 + B <=> BN + N 

N 2 + C <=> CN + N 

1.1 10 14 

46870 

23200 

22986 

4 

BN + B <=> B 2 + N 

CN + C <=> C 2 + N 

5.0 10 13 

31919 

13000 

19077 

5 

B 2 + N 2 <=> 2BN 

C 2 + N 2 <=> 2CN 

7.1 10 13 

14950 

5330 

3993 


Table 2. Enthalpy of Formation 14 


species 

A Hf, kJ/mol 

species 

A Hf, kJ/mol 

species 

A Hf, kJ/mol 

B 

560 ± 12 

c 

716.67 ± .46 

N 

472.68 ±0.1 

BN 

477 ± 125 

CN 

435.1 ± 10 



b 2 

829.7 ±33.5 

C 2 

837.7 ±3.8 

n 2 

0 


C. Material Response Model 

Two material response models have been tested. The first one simply allows the user to specify a 
uniform droplet surface temperature and a uniform vaporization mass flow rate as determined from 
experimental observation. The second one enables the computation of a spatially varying droplet 
temperature and mass flow from the droplet using a steady-state vaporization/melt approximation. 
An equilibrium vapor pressure is assumed for gaseous boron B(g) over the molten boron droplet 
B(l) at the tip of the solid boron fiber bundle. The boiling point of boron at p = lbar = 10 5 N/m 2 is 


7 of 23 


American Institute of Aeronautics and Astronautics 



T vap = 4137.895 K and the corresponding enthalpy of vaporization is A vap Hg = 480.344 kJ/rnol. 14 
(Uncertainties were not provided for these calculated values. For reference, the uncertainty in the 
heat of formation of boron gas is 12 kJ/rnol.) The Clausius-Clapeyron Equation is used to compute 
the equilibrium vapor pressure of boron as a function of droplet surface temperature. 


, PB , 2 
In 

PB, 1 


AvapHg ( 1 1 \ 

R \T2~tJ 


( 4 ) 


Table 3 provides a convenient reference for the boiling point of boron computed with Eq. 4 as a 
function of pressure levels used in this study. Chamber pressure in units of psig are used throughout, 
consistent with the experimental apparatus. Table 3 also shows these pressures in units of Pa and 
atmospheres. 


Table 3. Boiling Point of Boron 


p, psig 

p, Pa 

p, atm 

T, K 

0 

1.01325 10 5 

1 

4138. 

50 

4.46062 10 5 

4.4 

4629. 

100 

7.90800 10 5 

7.8 

4852. 

200 

14.8028 10 5 

14.6 

5122. 

400 

28.5928 10 5 

28.2 

5439. 

800 

56.1713 10 5 

55.4 

5808. 


A rigorous simulation would track the energy flux from the gas through the droplet and into the 
fiber bundle. A simplifying steady-state vaporization approximation is made to facilitate parametric 
studies. In the steady state vaporization approximation an energy balance across the liquid - gas 
interface of the droplet is evaluated. With reference to Fig. 4 the absorbed energy entering the 



Figure 4. A schematic of energy flux crossing the interface of the droplet surface and the atmosphere 
at the base of the plume. 


interface from above (the gas side) is the laser radiation aq ra d . The energy leaving the interface 
through the gas phase by molecular transport processes is q CO nv This term includes conduction 
and diffusion of multiple species in the gas phase. The energy re-radiated from droplet surface 
into free space is euT 4 . The energy carried by vaporization into the base of the rising plume is 
m g h w . The energy approaching the interface from the liquid below is carried by the melt mass 
flux rhmeithmeit ■ The energy leaving the interface in to the droplet is represented by q CO nd ■ In 
the steady state approximation we assume that the melt rate is balanced by the vaporization rate 
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so that the size of the droplet remains constant and there is no accumulation of energy in the 
droplet. Also, the conduction of energy into the droplet q CO nd is assumed to balance the melt flux 
of energy approaching the interface ihmelthmelt- All of these balances, represented by Eq. 5 with 
Qcond, canceling rh m eith me it, implicitly assume a one-dimensional variation of properties . 

aq ra d = + q conv + m g h w (5) 

In reality there are losses through the sides of the droplet that are not simulated in the idealized 
domain. Consequently, we will lump the effect of these losses into the parameter a which already 
includes the absorption of radiation from the laser at a non-orthogonal angle to the droplet surface. 
We ultimately expect that a is a parameter that can be adjusted to match observed recession rates 
of the boron fiber bundle. We must also confirm that the assumption of a constant droplet size 
(■ rhmeit = hig) can be achieved in the PVC process. 

V. Numerical Method 

Program LAURA 11,15 simulates the coupled radiation and ablation in the shock layer over 
a vehicle entering a planetary atmosphere at hypersonic velocities. This code was chosen as a 
starting point for simulation of the production rig because the algorithms for dealing with high 
temperature, reacting gas mixtures and tightly coupled ablation and radiation for carbon based 
ablators are already well established. It uses Roe’s averaging 16 with Harten’s entropy fix 17 and Yee’s 
Symmetric Total Variation Diminishing 18 algorithms for solving inviscid flux. Central differences 
are used for viscous flux terms. 

Simulation of the production rig required new modifications to the code including: (1) addition 
of boron species (B, BN, B 2 ) to thermodynamic and transport property data sets; (2) resetting a 
surface energy balance equation 19 under a laser radiation source based on equilibration of atomic 
boron vapor pressure with a liquid boron source; (3) adding buoyancy terms (gravitational field) to 
the momentum equations; and, (4) adding a porous wall boundary condition to mimic a pressure 
relief valve in the system to maintain constant pressure as mass and energy are added to the system. 
Items 1, 2, and 4 have already been discussed in the previous section on modeling. Item 3 is easily 
accommodated by adding the source term pg to the z-momentum equation. 

Conversion of a dependent variable set from species densities and energy to species partial 
pressures and energy have addressed some preconditioning issues in the application of a hypersonics 
flow solver to a nearly incompressible flow domain though some stability issues persist. These issues 
are addressed by running simulations with first-order time accuracy with At = 1. 10 -6 s. Ten sub- 
iterations are engaged before advancing to the next time step. Each sub-iteration uses implicit line 
relaxation extending from the lower wall and droplet surface to the upper wall. 

VI. Plume Simulations 


A. Parameter Space 

Nominal conditions for simulated PVC process are po = 200 psig, To = 300 K, po = 16.617 kg/nr 3 , 
in a molecular nitrogen atmosphere. As noted by Smith et al., 2 an elevated pressure in the PVC 
method is found to yield exceptionally long nanotube lengths (estimated to at least 1 mm based 
on ability to spin fibrils into 1 mm diameter yarn). Nanotubes may also be created at atmospheric 
pressure levels 20 but reported tube lengths at this condition are much shorter at 100 nrn. 

Reasons for longer nanotube length under pressure are not rigorously understood. Modeling 
of the nucleation process with molecular dynamics 21 and density functional theory 22 are required 
to understand the key environmental factors at elevated pressures. Wang and Yang 23 suggest 
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recently developed techniques of “nanothermodynamics” can also provide insight into the nucleation 
process. In all of these approaches the nucleation modeling must be informed by the macroscopic 
environment (temperature, pressure, species, residence time) provided here. 

Simulations of the macroscopic environment provide insights into the possible mechanisms 
through which the nucleation rates are enhanced. Certainly, the elevated pressures increase the 
collision rates of the component particles with the nucleation site. As will be seen in the subse- 
quent sections, raising the pressure levels also raises the level of BN molecules directly available 
for nucleation without the need to catalyze dissociation of N 2 in a liquid boron nucleation site. If 
condensation of BN molecules at the tip of a developing nanotube is a possible growth mechanism 
then it may prove to be a critical process to enable growth of exceptionally long BN nanotubes. 

The parameter space for which simulations are executed include chamber pressure levels of 0, 
50, 100, 200, 400, and 800 psig. The simulated chamber atmospheres include xn 2 = 1, Xn 2 = 0-5 
and XHe = 0.5, and xn 2 = 0.5 and XAr = 0.5. The power levels delivered to a 1 mm radius spot 
after subtracting losses include .05, 0.1, 0.2, and 2.0 kW. All of these simulations are intended 
to provide insight regarding their effect on the plume environment. The flow rates of BN gas in 
the plume are highlighted because of their potential importance to a tip growth mechanism. Only 
simulations using the steady state vaporization boundary condition are provided here in which mass 
loss rate and droplet temperature are computed as part of the boundary condition. 

B. Grid Convergence for Nominal Conditions 

The grid generation algorithm for the nominal condition simulations was described previously in 
the section Computational Domain, Grid, and Boundary Conditions. A grid convergence study is 
presented here to evaluate any effects of grid resolution on computed results. The fine grid and 
coarse grids use minimum mesh sizes that are a factor of two finer and coarser, respectively, than 
the nominal grid. The fine and coarse grids use a factor of two smaller and greater, respectively, 
time step of At = 1. 10~ 6 s used on the nominal grid simulations. The maximum allowed grid 
size, 2 mm, is unchanged among all three grids. Consequently, grid resolution far from the plume 
is equivalent. The study focuses on resolution at the plume base and across the radial edge of the 
rising plume. 



1 mm 


Figure 5. Contour lines of BN flow rates in the plume over an irradiated droplet of Boron are shown. 
Left side shows results for nominal grid. Right side shows results for fine grid. 
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Fig. 5 presents a composite contour plot of the nominal grid solution on the left and the fine 
grid solution on the right. The plume axis rises vertically in the center of the figure. The contours 
indicate the vertical flow rate of the BN species in units of mg/mm 2 -hr. (Nanotube production 
rates are most often quoted in units of mg/hr in the literature. Because the base radius of the 
plume base is 1 mm these units provide a convenient reference.) Note the concentric blue oval 
contour lines at the base of the plume. Negative values here indicate a recirculating flow in which 
the velocity is approaching the droplet surface. The foot of the plume is fed by cold, molecular 
nitrogen. As it encounters the thermal layer above the irradiated spot the temperature increases 
sufficiently to enable some dissociation of molecular nitrogen that reacts with atomic boron rising 
from the droplet. A slight, unsteady pulsing of the plume is evident early in its evolution on the 
fine grid. The pulsing is not evident in the nominal grid simulation and eventually damps out on 
the fine grid. (Unsteady motion persists at higher absorbed laser power levels as will be discussed 
subsequently.) 

Fig. 6 presents the centerline variation of plume properties as a function of distance above the 
droplet. A log scale is used to focus on chemistry at the plume base. The coarse grid solution 
shows higher temperature and BN mass fraction at the base of the plume compared to the nominal 
and fine grid solutions. The coarse grid solution clearly dissipates the cbn and temperature profile 
as the plume rises. In contrast, grid convergence is substantially achieved at the base of the plume 
for z < 0.1 mm across the recirculating flow as judged by agreement of the medium and fine grid 
results. Fig. 6(d) shows the evolution of the temperature profile at 5 ms intervals from the fine 
grid. The blue line in this figure shows a profile very early in the development of the plume at 0.2 
ms. The red line in this figure indicates the medium grid solution at 40 ms. For 0.1 < z < 1 mm 
the medium grid solution starts to dissipate the higher curvature elements of the fine grid profiles. 
For z > 1 mm the plume growth on the medium grid starts to lag the growth on the fine grid as 
evident by the higher temperature profiles for a given time. 

All observations relating to temporal evolution of the plume should be considered qualitative 
rather than quantitative. First, the axi-symmetric domain and boundary conditions restrict tem- 
poral evolution - the core of the plume is constrained to the centerline of the domain. Second, the 
assumption of a steady-state vaporization boundary condition requires a finite amount of time to set 
up. Given these caveats one notes that contour levels in the base of the plume set up very quickly. 
Consider, for example the line plot evolution of temperature in Fig. 6(d) where each line indicates 
the extent of the temperature profile in 5 ms increments. After only 15 ms the temperature at z = 
1 mm is within 5% of its value at 40 ms. Species mass fractions follow the same trend. 

A check of cell Reynolds numbers Re ce u = pwAx/p across the diffusion layer of the plume 
on the nominal grid indicate peak values as a function of pressure of 1.1 at 0 psig, 38 at 50 psig, 
62 at 100 psig, 94 at 200 psig, 137 at 400 psig, and 241 at 800 psig. Given the small differences 
between the nominal and fine grid solutions at 200 psig for z > 0.1 mm it is assumed that numerical 
dissipation affects the solution for at least Re ce u > 90. For z < 0.1 mm the diffusion layer is better 
resolved with at least an order of magnitude reduction in cell size at the foot of plume where gas 
flow is predominantly in the radial direction. Evidence of grid convergence in the plume base is 
critical because the conditions here drive subsequent plume evolution and define maximum available 
reservoir levels of BN. 

C. Chamber Pressure Effects 

Comparisons of the plume shape and properties at 0 psig (atmospheric pressure) and 200 psig (14.6 
atmospheres) are presented in Figs. 7 and 8. The plume radius expands from the base at 0 psig 
and contracts at 200 psig for the same delivered power level of .1 kW. The vertical black line in 
these figures to the left and right of the axis indicate the 1 mm radius corresponding to the radius 
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(a) mass fraction of BN, t = 0.04 s (b) flow rate of BN, t = 0.04 s 


At =5. 10 7 S 




(c) temperature, t = 0.04 s (d) temperature evolution on finest grid 

Figure 6. Centerline variation of properties in a rising plume of gas over an irradiated droplet of 
Boron as function of grid and time. 
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of the irradiated spot on the boron droplet on the floor of the domain. Peak temperature at the 
base and peak mass fraction of BN are higher for the 200 psig case compared to the 0 psig case. 
However, the BN mass fraction at 200 psig declines more quickly. The corresponding flow rates 
of boron derived species rising to the condenser loop (not shown) are presented for BN in Fig. 
8(a) and for the sum of B and B 2 in Fig. 8(b). There is no recirculation pattern at the base of 
the plume at 0 psig. Results indicate that BN returns to B 2 and N 2 as the plume rises and cools 
under assumption of equilibrium chemistry and under the assumption of single phase flow. The 
single phase flow assumption forces a supersaturated level of BN and B in the rising plume. The 
equilibrium flow assumption may deplete BN faster than collisional processes would allow because 
the availability of suitable atomic collision partners is diminished in the rising plume as it cools. 
This question will be revisited in the next section on Gas Chemistry Model Effects. 
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Figure 7. Comparison of plumes at 0 psig (left) and 200 psig (right). 




(b) Flow rate of B and B 2 


Figure 8. Comparison of Boron vertical flow rates at 0 psig (left) and 200 psig (right). 


Figures 9-11 are intended to capture maximum values of key plume environmental variables 
that tend to drive plume dynamics. These maximum values occur at or near the base of the 
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plume where grid resolution is thought to be adequate (based on the earlier Grid Convergence 
discussion) to establish the trends associated with chamber pressure. The figures show results as 
a function of total chamber pressure on the left and as a function of partial pressure of nitrogen in 
the chamber on the right. Helium and argon are added to the parametric study here to investigate 
if the presence of a lighter or heavier noble gas may influence the environment in a way to enhance 
the flow rate of boron-richgases. (The results shed no light regarding how the presence of noble 
gases may effect nucleation rates.) Isolating the effect of nitrogen partial pressure is intended to 
capture the importance of the availability of nitrogen as a reaction partner with boron at the base 
of the plume. There is no evident advantage to diluting the nitrogen atmosphere with a noble gas 
based on these results. 

The figures show a clearly defined peak for BN mass fraction and flow rate between 50 and 200 
psig. The maximum values are highest for a pure nitrogen atmosphere though the trends are evident 
in all of the tested atmospheres. It is thought that rapidly rising temperature between 0 and 100 
psig accounts for a greater availability of atomic nitrogen which dissociates at temperatures above 
4000 K. Beyond 200 psig the increase in peak temperature tapers off (Fig. 11) while dissociation of 
nitrogen is suppressed for a given temperature as the pressure rises. Consequently, the production 
of BN diminishes beyond 200 psig. If nucleation from a supersaturated flow of BN is a significant 
driver to the BNNT production rate or nanotube length then this local maximum is an important 
process design environmental variable. 




(a) total chamber pressure 


(b) nitrogen partial pressure in chamber 


Figure 9. Variation of maximum BN mass fraction as a function of pressure and atmospheric con- 
situents. 


D. Gas Chemistry Model Effects 

Earlier sections on Gas Properties and Chamber Pressure Effects discussed the concern that an 
equilibrium assumption may cause BN gas to be converted to B 2 and N 2 at an unnaturally high 
rate because of the limited availability of atomic collision partners. (This issue is related to the 
need for catalytic converters to break down exhaust gases from an automobile.) Two chemical 
kinetic models for the B-N system are proposed based on the C-N system as documented in Table 
1. Kinetic model 1 includes all 5 reactions listed in Table 1. Kinetic model 2 includes only the 
first 4 reactions listed in Table 1. It omits the final reverse reaction, effectively a double shuffle 
reaction, in which two molecular boron nitride molecules collide, dissociate, and then recombine 
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Figure 10. Variation of maximum vertical flow rate 
consituents. 



(b) nitrogen partial pressure in chamber 

of BN as a function of pressure and atmospheric 


into diatomic boron and diatomic nitrogen. Without this reaction, the only path to dissociate 
BN requires an atomic intermediary as either a collision partner (reactions 3 and 4) or a product 
(reaction 2). The template for reaction 5 had originally been proposed to explain dissociation of 
CN in a shock tube 24 and later to explain measurements behind reflected shock waves 25 but no 
such mechanism is used in kinetic models for air (O-N) systems. 

Figure 12 compares mass fraction, BN flow rate and integrated flow rates across the plume as 
a function of height. In Fig. 12(a) the two kinetic models show identical levels of mass fraction 
cbn i n the lower portion of the plume (z < 1 mm) extending through the recirculating flow where 
BN is initially formed. The equilibrium model predicts approximately 10% lower values of cbn in 
this same region. For z > 1 mm kinetic model 1 approaches the equilibrium model while kinetic 
model 2 shows that BN persists as the plume continues to rise and cool. This trend is also evident 
in the vertical flow rate of BN in Fig. 12(b) where significant levels of BN persist beyond 2 cm to 
the top of the plume. The flow rate per unit area is integrated across the area of the plume out 
to 1 mm radius to show the flow rate of BN in mg/hr in (c). An equivalent integration produces 
the combined flow rate of B and B 2 in mg/hr. The BN flow levels using kinetic model 2 are small 
compared to B and B 2 but the rate is roughly equivalent to estimated production rates of BNNTs 
in the PVC process. Kinetic model 2 results have no experimental validation. Its formulation using 
the C-N system template is thought to provide a credible example of how an actual flow would 
behave. 

E. Laser Power Effects 

The effects of absorbed laser power to a 1 mm radius spot on plume development are presented in 
Figs. 13 and 14. Simulated absorbed power levels include .05 kW, 0.1 kW, 0.2 kW, and 2.0 kW in 
a 200 psig chamber with \n 2 = 1- Kinetic model 2 is employed in all simulations. The profiles are 
captured at 80 ms. 

As absorbed laser power increases the temperature at the base of the plume (Fig. 13(a)) ap- 
proaches the boiling point of boron as documented in Table 3. The boiling point is achieved at 
2.0 kW and the blowing rate of boron completely displaces nitrogen at the bottom of the plume. 
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(a) total chamber pressure (b) nitrogen partial pressure in chamber 

Figure 11. Variation of maximum temperature as a function of pressure and atmospheric consituents. 


The sum of eg and cg 2 equals 1 for 2.0 kW absorbed power (Figs. 13(c) and (d)) and cbn goes to 
zero (Figs. 13(b)) because no nitrogen is available to react with boron at the droplet surface. (The 
boiling droplet would not stay suspended from the fiber bundle in the current configuration - it is 
presented here as an interesting limiting condition.) A dip in the temperature profile at 2 = 0.1 
mm for the two lowest power levels is indicative of a recirculation pattern at the base of the plume. 
As power is increased, the blowing rate increases sufficiently to displace the recirculating flow. 

Profiles of the flow rates of boron species are presented in Fig. 14 for each of the four power 
levels. The centerline rates (Fig.l4(a-b)) show the convective flux of boron containing species as a 
function of height. The integrated flow rate (Fig.l4(c-d))plots the integral of the rising convective 
flux of boron containing species from x = 0 (centerline) to x = 1 mm capturing the vertically 
convecting mass flux in the entire plume as a function of height. Note that these plots include 
only the vertical component of convective flux. This component represents almost 100% of the 
boron flux in the upper portion of the plume. If the plume is unsteady, there can be a significant 
horizontal component of flux. More importantly, for low blowing rates at lower absorbed power 
levels, the diffusive component of boron mass flux exceeds the convective component at the base of 
the plume. 

The dashed lines in Fig. 14(d) denote the vaporization rate of boron from the droplet surface. 
They are extended across the figure so that one may compare the instantaneous mass loss rate 
at 0.08 s to the vertical mass flux rate in the established plume. For the two lowest laser power 
cases agreement is excellent for z > 1 mm. For the two highest power cases the simulated plume is 
unsteady, consequently the integrated flow rate shows variation as a function of height caused by 
pulsing as a function of time. The difference between the integrated convective flow rate and the 
vaporization rate for z < 0.01mm indicates the level of diffusive flux of boron leaving the droplet. 

The sustained integrated flow rate of BN in the upper part of the plume is approximately 3% 
of the vaporization rate for all of the power levels at 200 psig. These sustained levels essentially 
disappear if BN can form B 2 and N 2 at low temperatures without atomic intermediaries as forced 
by the chemical equilibrium approximation or as allowed in a kinetic model including reaction 5 of 
Table 1. These sustained levels could rise or fall depending on the kinetics of the real B-N system 
compared to the proposed system simulated here. 

The computed integrated vaporization rate fh for a feeder fiber bundle with density p ro d and 
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Kinetic Model 1 




(a) mass fraction BN (b) vertical flow rate of BN 




(c) integrated flow of BN across plume 


(d) integrated flow of B and B 2 across plume 


Figure 12. Centerline variation of properties in a rising plume of gas over an irradiated droplet of 
Boron as function of gas chemistry model. 
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T, K 




(a) temperature (b) mass fraction of BN 




(c) mass fraction of B (d) mass fraction of B 2 

Figure 13. Centerline variation of thermodynamic properties in a rising plume of gas over an irradiated 
droplet of Boron as function of laser power. 
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(a) centerline flow rate of BN (b) centerline flow rate of B and B 2 




(c) integrated flow rate of BN across plume 


(d) integrated flow rate of B and B 2 across plume 


Figure 14. Variation of centerline and integrated production rates in a rising plume of gas over an 
irradiated droplet of Boron as function of laser power. Dashed lines indicate blowing rate at surface. 
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radius r ro d is easily transformed to the recession rate of the feeder fiber bundle z ro d with Eq. 6. 

fh' — ZrodProd'KT'rod ( 6 ) 

Assuming a boron fiber bundle density p ro d of 2.46 g/cm 3 the vaporization rates in Fig. 14(d) 
transform to recession rates of 1.33 cm/hr for 0.05 kW, 13.5 cm/hr for 0.1 kW, 93.8 cm/hr for 
0.2 kW and 1393 cm/hr for 2.0 kW. For reference, the corresponding heating rates vary from 1592 
W/crn 2 to 63662 W/crn 2 from the lowest to highest absorbed laser power levels. The fiber bundle 
recession rate will serve as a relatively simple metric to calibrate absorbed power after losses. Its 
utility will depend on ability to achieve a steady state vaporization condition without mechanical 
loss of droplets from the end of the feeder fiber bundle. 

VII. Experimental Observations versus Simulations 

PVC process development has only recently included the simulation component described 
herein. Furthermore, the ability to execute these simulations has become available during a down 
time in which the production rig is being moved from the Department of Energy, Jefferson Lab in 
Newport News, Virginia to NASA Langley Research Center in Hampton, Virginia. Consequently, 
there has not yet been a concerted effort to gather experimental measurements for validation of the 
simulations. The BNNT apparatus became operational again at NASA Langley on May 17, 2012 
and such validation experiments are currently in the planning stages. Two key production process 
observations 26 are noted and discussed here. 

It is observed that significant condensate forms on the condenser loop at locations from approx- 
imately 1 cm to 7 cm above the droplet surface. Nominal grid solutions indicate that it takes 0.16 
seconds for the top of the plume to rise to a height of approximately 5 cm. The plume is still slowly 
rising at this point in time. Plume simulations show greater attained height at equivalent times 
with greater delivered laser power (lower assumed losses). However, with the current modeling of 
steady state ablation from a spot on the floor of the chamber these higher power simulations indi- 
cate a greater recession rate than estimated from observations. Small perturbations to the current 
models may make the plume rise faster - for example, reducing the estimated laser spot diameter 
while increasing delivered power. A better solution is to increase fidelity of the laser - droplet - fiber 
bundle interaction region by expanding the simulation to include the domain around the suspended 
droplet. This modification exposes more droplet surface to feed the plume and allows an upwelling 
of nitrogen from below the droplet to potentially assist loft. 

It is observed that the “cotton-like” condensate produced at a rate of 50 to 100 mg/hr contains 
up to 50% nanotubes by mass. The remaining condensate includes BN particle flakes and solidified 
boron droplets. Nanotubes are most often seen extending from these droplets suggesting that liquid 
boron droplet sites are critical to the nucleation of BNNT nanotubes. In the process described by 
Arenal et al. 27 and Lee et al. 20 it is assumed that liquid boron assists in the catalysis of molecular 
nitrogen to build the BN nanotube from the root. The simulations here suggest that a supersatu- 
rated state of BN gas in the plume provides a reservoir from which extremely rapid condensation of 
nanotubes progresses when a “suitable” nucleation site is encountered. The saturated equilibrium 
vapor pressure of boron with liquid boron, ps, S at, and the saturated equilibrium vapor pressure of 
BN with hexagonal boron nitride solid, p bn, sat, are calculated as a function of temperature using 
Eq. 4. Estimated reference triple point conditions used in Eq. 4 for BN are pbn = 400 ± 20 
Pa, 28 Ten = 3400 ± 20 K, 28 and Hg N sub = 669.44 kJ/rnol. 14 The reference conditions for B were 
already discussed in the section on Material Response Model. The local vapor pressure of B and 
BN non-dinrensionalized by the corresponding saturated vapor pressure is presented in Fig. 15. 
The supersaturation factors exceed 100 for z > 1 mm. The simulated mass flux levels of BN, B, 
and B 2 are sufficiently high to produce the observed levels of condensate formed in the chamber. 
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z, m 

Figure 15. Profiles of B and BN vapor pressure non-dimensionalized by the corresponding saturation 
vapor pressure as a function of distance above the droplet at po = 200 psig using kinetic model 2. 


Finally, we note that as chamber pressure is increased above 100 psig for absorbed laser power 
equal to 0.1 kW a recirculation region forms at the base of the plume. Recirculation provides 
extended residence time in a BN-rich environment for any particles (potential nucleation sites) 
that enter this region. As noted previously, higher absorbed power levels tend to increase the 
vaporization rate and displace any recirculation pattern. 

VIII. Summary 

Flow in a pressurized, vapor condensation (PVC) boron nitride nanotube (BNNT) production 
rig is modeled within Program LAURA. A laser provides a thermal energy source to the tip of 
a boron fiber bundle in a high pressure nitrogen chamber. The thermal energy source causes a 
molten droplet to form at the end of the fiber bundle. The boiling point temperature of boron 
increases with an increase in chamber pressure. The droplet temperature also rises with increasing 
chamber pressure for a given absorbed power level under assumptions of equilibrium vapor pressure 
and steady state vaporization. A plume of boron-rich gas rises from the droplet surface. The 
buoyancy driven flow is modeled as mixture of thermally perfect gases (B, B 2 , N, N 2 , BN) in 
either thermochemical equilibrium or chemical nonequilibrium assuming steady-state melt and 
vaporization. 

An idealized, axisymmetric chamber is used to model the PVC process. The droplet is positioned 
on the chamber floor at the axis, a 1 mm radius spot is irradiated, and a time-dependent evolution 
of the plume is simulated. Parametric simulations involving chamber pressure, equilibrium and 
non-equilibrium gas models, laser power levels, and dilution of nitrogen with noble gases in the 
chamber are executed. The simulations are intended to define the macroscopic thermochemical 
environment from which boron-rich species, including nanotubes, condense out of the plume. Only 
the gas phase is considered. Flow rates of BN vapor are monitored on the assumption that they 
could contribute to enhanced nucleation rates of BNNTs. The simulations for a nominal condition 
defined by po = 200 psig and absorbed laser power of 0.1 kW to a 1 mm radius spot on a molten 
boron droplet indicate: 

1. A high temperature environment (T > 4400K) forms within 1 mm of the surface at elevated 
chamber pressures (po > 50 psig) sufficient to dissociate molecular nitrogen and form BN at 
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the base of the plume. 

2. BN vapor flow rate exhibits a maximum at the plume base for po somewhere between 50 and 
200 psig. 

3. Under an equilibrium chemistry model BN is transformed to B 2 and N 2 as gas rises and the 
plume cools. 

4. Under a chemical kinetic model BN vapor flow levels are frozen as the plume continues to 
rise. The BN flow rate is approximately 3% of the net vaporization rate at po = 200 psig for 
absorbed laser power varying from 0.05 kW to 2.0 kW 

5. A recirculation region begins to form at the base of the plume at po = 100 psig. Recirculation 
provides any existing particulates an opportunity for extended residence time in a BN rich 
environment, even with the equilibrium gas model. 

The simulations indicate that the high pressure environment causes significant levels of BN 
vapor in the plume that are not present at atmospheric pressure. We believe that a supersaturated 
reservoir of BN vapor is a significant contributor to the formation of high quality BNNTs produced 
in the PVC process. This theory still requires experimental validation that the physical models 
used in the simulation of the macroscopic environment of the plume are sufficiently accurate. The 
theory also requires simulations of the nucleation process itself to evaluate conditions that favor 
formation of nanotubes as opposed to the formation of amorphous particles of boron and boron 
nitride. 
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